% 1) construct m = {m(tj)}, a matrix where the rows are funds and the
% columns are the SDF corresponding to each of the J cash flows of the fund
% (note: these are padded with zeros to make a balanced matrix)
% 2) compute gt = m*r', where r={r(tj)} i.e. matrix with fund cash flows
% 3) compute dg = 

% r = t x nq x J
% x = t x nx x J 

function [gt, dg] = fundMCexpaff_hvar(b,v,x,r,p,h,fast)   %h is only there to be consistent with other functions that need h

t = size(r,1);      % # funds
nq = size(r,2);     % number of observed assets to be priced (e.g., fund mkt rf) 
nx = size(x,2);     % # pricing factors (e.g. const mkt) 

J = size(x,3);      % max # cash flows in fund 

mm = zeros(t,1,J);

x(isinf(x)) = 0;    % these are the periods post end-of-fund that were padded with zeros. 
x(isnan(x)) = 0;    % ... and VFfund_CF always zero for those obs

if v == -1          % no factor (co)variance matrix specified 
    for j = 1:J         % # cycle through all fund cash flow periods   
        mm(:,1,j) = exp( squeeze(x(:,:,j)) * b );      % sdf t x 1 x J
    end
else                % use factor covariance for each horizon. Horizons presumed to be first factor, i.e., x(:,1,:), and log risk-free rate the second factor.
    if nx == 3      % if only one factor (besides horizon and Rf)
        bb = [zeros(t,1) ones(t,1)*(b-1) ones(t,1)*(-b)];  % first entry has horizon zero so first parameter does not matter
        for j = 1:J         % # cycle through all fund cash flow periods   
            vsq = squeeze(v);   % since v is 1x1x15
            if j > 1
                vind = min(15, ceil(squeeze(x(:,1,j))) );   % index to correct horizon for each cash flow
                vind(vind==0) = 1;                          % if horizon is exactly zero then this doesn't matter (just to avoid error messages)
                bb = [vsq(vind)*b*(b-1)/2 ones(t,1)*(b-1) ones(t,1)*(-b)];
            end
            mm(:,1,j) = exp( sum(squeeze(x(:,:,j)) .* bb, 2) );      % sdf t x 1 x J    
        end
    elseif nx == 4      % if two factors (besides horizon and Rf)        
        
        bb = [zeros(t,1) ones(t,1)*(b(1)-1) ones(t,1)*(-b(1)) ones(t,1)*(-b(2))];                       
        
        
        % precompute possible first column entries for parameters (for speed)
        bb1 = zeros(15,1);
        for k = 1:15
            vsq = squeeze(v(:,:,k));
            bb1(k) = -b'*diag(vsq)/2 + b'*vsq*b/2;
        end        
        % # cycle through all fund cash flow periods   
        for j = 1:J         
            if j > 1
                vind = min(15, ceil(squeeze(x(:,1,j))) );   % index to correct horizon for each cash flow
                vind(vind==0) = 1;                          % if horizon is exactly zero then this doesn't matter (just to avoid error messages)                                            
                bb(:,1) = bb1(vind);
            end
            mm(:,1,j) = exp( sum(squeeze(x(:,:,j)) .* bb, 2) );      % sdf t x 1 x J    
        end               
    else
        disp('Can''t do more than 2 factors!');
    end
end
m = repmat(mm,[1 nq 1]); %to have same dim as cash flow array r: t x nq x J
gt = sum(bsxfun(@times,r,m),3) - p;    % t x nq

 
dg = zeros(nx,nq); 

if fast == 0 || nargin < 6
   dgt = zeros(nx,nq);  
   for i = 1:t 
      if nx > 1  
         adddg = squeeze(x(i,:,:))*(squeeze(m(i,:,:)).*squeeze(r(i,:,:)))';
      else 
         adddg = permute(x(i,:,:),[1 3 2])*(squeeze(m(i,:,:)).*squeeze(r(i,:,:)))'; 
      end
      dgt = dgt+adddg;
   end
   
   dg = dgt/t;

end

%just for debugging
%global sdf
%sdf = squeeze(mm); 
%global gs
%gs = gt; 

%{
size(dg)
size(cov(gt))
W = eye(3);
inv(dg*dg')
e = eig(dg*W*dg');
e(1)
e(2) 
dg*W*cov(gt)*W*dg'
W(3,3) = 0.00001;
dg*W*dg'
e= eig(dg*W*dg');
e(1)
e(2) 
ddg = dg(:,1:2); 
e= eig(ddg*ddg');
S = cov(gt);
e(1)
e(2)
(1/t)*inv(ddg*ddg')*ddg*S(1:2,1:2)*ddg'*inv(ddg*ddg') 
%}
